Astronomy & Astrophysics manuscript no. ComprSGSModel 


©ESO 2010 


October 22, 2010 





A Fluid-Dynamical Subgrid Scale Model for 
Highly Compressible Astrophysical Turbulence 

W. Schmidt 1 and C. Federrath 2 3 - 4 

1 Institut fur Astrophysik, Universitat Gottingen, Friedrich-Hund-Platz 1, D-37077 Gottingen, Germany 
e-mail: schmidt@ast.ro . physik . uni -goettingen . de 

2 Zentrum fur Astronomie der Universitat Heidelberg, Institut fur Theoretische Astrophysik, Albert-Ueberle-Str. 2, D-69120 
Heidelberg, Germany 

e-mail: chf eder@ita . uni-heidelberg . de 

3 Max-Planck-Institut fiir Astronomie, Konigstuhl 17, D-691 17 Heidelberg, Germany 

4 Ecole Normale Superieure de Lyon, CRAL, 69364 Lyon Cedex 07, France 

ABSTRACT 

Context. Compressible turbulence influences the dynamics of the interstellar and the intergalactic medium over a vast range of length 
scales. In numerical simulations, phenomenological subgrid scale (SGS) models are used to describe particular physical processes 
below the grid scale. In most cases, these models do not cover fluid-dynamical interactions between resolved and unresolved scales, 
or the employed SGS model is not applicable to turbulence in the highly compressible regime. 

Aims. We formulate and implement the Euler equations with SGS dynamics and provide numerical tests of an SGS turbulence energy 
model that predicts the turbulent pressure of unresolved velocity fluctuations and the rate of dissipation for highly compressible 
turbulence. 

Methods. We test closures for the turbulence energy cascade by filtering data from high-resolution simulations of forced isothermal 
and adiabatic turbulence. Optimal properties and an excellent correlation are found for a linear combination of the eddy-viscosity 
closure that is employed in LES of weakly compressible turbulence and a term that is non-linear in the Jacobian matrix of the 
velocity. Using this mixed closure, the SGS turbulence energy model is validated in LES of turbulence with stochastic forcing. 
Results. It is found that the SGS model satisfies several important requirements: 1. The mean SGS turbulence energy follows a power 
law for varying grid scale. 2. The root mean square (RMS) Mach number of the unresolved velocity fluctuations is proportional to the 
RMS Mach number of the resolved turbulence, independent of the forcing. 3. The rate of dissipation and the turbulence energy flux 
are constant. Moreover, we discuss difficulties with direct estimates of the turbulent pressure and the dissipation rate on the basis of 
resolved flow quantities that have recently been proposed. 

Conclusions. In combination with the energy injection by stellar feedback and other unresolved processes, the proposed SGS model 
is applicable to a variety of problems in computational astrophysics. Computing the SGS turbulence energy, the treatment of star 
formation and stellar feedback in galaxy simulations can be improved. Further, we expect that the turbulent pressure on the grid scale 
affects the stability of gas against gravitational collapse. The influence of small-scale turbulence on emission line broadening, e. g., of 
O VI, in the intergalactic medium is another potential application. 

Key words. Hydrodynamics - Turbulence - Methods: numerical - ISM: kinematics and dynamics - Galaxies: evolution 

1 . Introduction If the unresolved turbulent velocity fluctuations reach a non- 
negligible fraction of the speed of sound, they give rise to a tur- 
bulent pressure in addition to the thermal pressure of the gas. 

The effects of numerically unresolved turbulence have recently This contribution to the pressure is proportional to the energy 

met increasing atte ntion in a variety of astrophysica l simulations density of the SGS turbulence. Turbulent pressure effects in the 

H \ (see, for instance , Scannapi eco & Briiggenl l2008t iMaier et all b aryonic gas compone nt of star-forming galaxies are discussed 

l2009t|joung et al.l f2009; Oppenheimer & Dave 2009). Some ap- in lBurkert et all d2009l) . In contemporary numerical simulations 

proaches comprise subgrid scale (SGS) models, although these of disk galaxies (e g.. | Dobbs et al.ll2008t iTasker & Tanl f2009: 

models are basically phenomenological parameterizations of as- lAgertz et aT1l2009L 120101) or in galactic-scale simulations of the 

trophysical processes on length scales smaller than the grid inters tellar medium (e. g., liming & Ma c Lowll2006t [joung et all 

scale. The full multi-scale dynamics of turbulence, however, l2009l) . the minimal grid scale (or the SPH smoothing length) 

is not embraced. The essence of an SGS model in the fluid- A > 1 pc. Since this length scale is comparable to the size 

dynamical sense is that, at high Reynolds numbers, energy is of molecular clouds, the unresolved turbulent velocity fluctu- 

transported through a turbulent cascade from larger, numeri- ations can exceed the speed of sound in the cold gas phase, 

cally resolved length scales to the subgrid scales. The energy of Consequently, it can be expected that a significant turbulent pres- 

the unresolved turbulent velocity fluctuations is eventually dis- sure is caused by turbulence below the grid scale. To compute 

sipated into heat. Numerical simulations, in which an explicit the turbulent pressure, which has an impact on the star formation 

closure for the turbulence cascade is applied on the grid scale, rate through the stability of the gas against gravitational collapse, 

are called large eddy simulations (LES). A closure is an approx- a model for the highly compressible regime is indispensable, 
imation to an SGS quantity in terms of resolved flow quantities. 
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Bona zzola et all (Il987l) and iBonazzola et alJ (1 19921) formulate 
an analytic theory to calculate the turbulent pressure of isotropic 
compressible turbulence on the integral length scale. Applying 
an SGS model, on the other hand, the turbulent pressure can 
be computed for a ny length scale within the inertial subrange. 
IJoung et al] J2009) present an SGS model that is based on the 
equation for the kinetic energy of the unresolved turbulent ve- 
locity fluctuations, the so-called SGS turbulence energy, where 
energy is solely supplied by supernova feedback. Since the non- 
diagonal SGS turb ulence stresses are neglected, the model of 
IJoung et al] (12009) reduces the effects of SGS turbulence to the 
turbulent pressure alone, and the turbulence energy cascade, i. e., 
the production of SGS turbulence by the shear of the numerically 
resolved flow, is not considered. 

A furt her example for this type of SGS models is the 
model of Scanna pieco & Bragged (120081) for the simulation 
of Rayleigh-Taylor-driven turbulence in active galactic nu- 
clei, where it is assumed that SGS turbulence is produced by 
buoyancy processes on unresolved length scales only. These 
processes are modelled by an equation for the characteris- 
tic length scale of the Rayleigh-Taylor instability. In contrast, 
ISchmidt et al.l (l2006bl) incorporate unresolved buoyancy effects 
into an SGS model that includes the production by shear for the 
treatment of turbulent combustion in thermonuclear supernovae. 

In the cosmo logical simulations of galaxy clusters by 
Maier et aH d2009l) . the role of SGS turbulence has been explored 
with a numerical technique that combines adaptive mesh refine- 
m ent and LES. They app ly the SGS turbulence energy model 
of ISchmidt etall (l2006ah . The main effect of the SGS model 
is an enhancement of the turbulent heating in the cluster core. 
The SGS turbulence energy also serves a s a tracer of turbulenc e 
production in the intergalactic medium (Iapichinoet al . 2010). 
However, a deficiency in these simulations is that the employed 
SGS model is only applicable to moderately compressible turbu- 
lence. Shocks are treated tentatively, i. e., SGS turbulence pro- 
duction is suppressed in the vicinity of shock fronts. While this 
is not a severe constraint for the bulk of the intracluster medium, 
in which the Mach numbers of the turbulent flow are small com- 
pared to unity, an erroneous production of SGS turbulence en- 
ergy is likely to occur near accretion shocks in the outer regions 
of the cluster. 

In this article, we improve on the previous approaches to 
SGS modelling by addressing the closure problem for highly 
compressible turbulence. In Sect. [2 we discuss the meaning of 
the compressible Euler equations in the context of computational 
fluid dynamics. The verification of the proposed closure and the 
calibration of the closure coefficients are presented in Sect. [3] 
Data f rom several high-resolution simulations of forced turbu- 
lence (ISchmidt et al. 2003, l2009t lFederrathetaT]l2010bl) allow 
us to compute the rate at which energy is transferred from length 
scales greater than the filter length to smaller length scales and 
to test the correlation with different closures. As a result, we 
propose a combination of the eddy-viscosity closure, which has 
successfully been used in LES of incompressi ble turbulence, 
and a non-linear closure that is put forward bv I Woodward et ail 
(2006). Then we show that physically reasonable statistics of the 
SGS turbulence energy and the rate of dissipation are obtained 
for varying grid resolutions and forcing in LES of supersonic 
turbulence (Sect. [4]). Furthermore, we investigate correlations of 
the SGS quantities with quantities derived from the numerically 
resolved flow. We demonstrate that the turbulent pressure and 
energy dissipation cannot be predicted in a straight-forward way 
on the bas is of the resolv ed turbu lent flow, as pr oposed, for in- 
stance, bv lPan et alJ d2009l) and lZhu et alJ ((2010). Instead, a full 



SGS model is needed to estimate unresolved turbulence effects. 
In the last Section, we summarize the results and discuss po- 
tential astrophysical applications of the closure for the highly 
compressible turbulence cascade in combination with the phe- 
nomenological approaches described above. 

2. The compressible Euler equations with 
subgrid-scale dynamics 

The Reynolds number of turbulent flows in astrophysics is usu- 
ally considered to be high enough so that the approximation of 
an inviscid fluid can be applied on numerically resolvable length 
scales. For a physically complete picture, we begin with the 
compressible Navier-Stokes equations encompassing all physi- 
cal length scales. This acknowledges the fact that perfect fluids 
do not exist in nature and that the notion of viscous dissipation 
is essential for turbulence. The fluid-dynamical variables deter- 

oo 

mined by this set of equations are denoted by g for the mass den- 
sity of the gas, u for the velocity of the flow, etcQThe resolution 
of a numerical simulation is given by the size of the grid cells A, 
which is called the cutoff scale or the grid scale. A consistent for- 
mulation of the equations of fluid dynamics with a cutoff scale A 
can by derived from the Navier- Stokes equ a tions b y means of the 
filter formalism introduced by iGermanol {1992). Generalizing 
this f ormalism to compre ssible fluid dynamics is straightforward 
(see ISchmidt et al.ll2006al) . The basic idea is to identify the nu- 

co 

merically computed solution with filtered variables g := (p)a, 

coco 

u := (gu)k/g, etc. The filter operator (-)a smoothes the physi- 
cal variables that are given by the Navier-Stokes equations over 
the length scale A. In LES, the filtering corresponds to the dis- 
cretization of the equations of fluid dynamics. The dynamical 
equations for the computable, filtered quantities are similar to 
the unfiltered equations, with additional terms that are related to 
the subgrid-scale dynamics on length scales I < A. 

Let us consider the dynamical equation for the momentum 
density of the fluid, which is given by the partial differential 
equation (PDE) 



8 /coco\ _ /cooo co\ do /co °°\ _< x > 

- \ou\ +V \pu®uj = p\g + fj-VP 



+ V • o% 



(1) 



where g and / are the accelerations due to gravity and other me- 

co 

chanical forces acting on the fluid, and P is the thermal pressure. 

CO 

The viscous dissipation tensor <x is defined by 



^ CO 

<r u =2vQ\Sij--d6ij 



(2) 



where v is the microscopic viscosity of the fluic0, the rate of 

°° co oo 

strain Sij is the symmetic part of the Jacobian matrix u ul - djUt, 

CO 

and d = m,-,,. Applying a homogeneous filter operator that is uni- 
form in time, Eq. (Q]) is converted into an equation for the fil- 

coco 

tered momentum, gu — (gu). This equation has the same form 
as the original equation, except for one term. Because of the non- 
linear advection term, the filtering introduces a stress term that 



1 A mathematical proof for the existence of a solution is still on the 
agenda of the Clay Mathematics Institute Millennium Prize Problems, 
even in the case of incompressible turbulence. 

2 Although we consider compressible fluid dynamics, for brevity, we 
neglect the second viscosity that is related to the divergence of the flow. 
This does not affect subsequent arguments. 
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accounts for the interaction between the numerically resolved 
flow and velocity fluctuations on the subgrid scales: 

—(gu) + V • (gu ® u) = p(g + /) - VP + V • (<r + T sgs ) , (3) 
where the SGS turbulence stress tensor is defined by 

OOOO CO 

i"sgs = -{q u ® u )a + gu®u. (4) 
In the following, the components of r sgs are simply denoted by 

coco oo 

Tij. The second-order moment {gu ® h)a is not explicitly com- 
putable in LES because the variations of the mass density g and 

CO 

the velocity u below the grid scale are unknown. For this rea- 
son, an approximation in terms of filtered quantities has to be 
devised. This is the closure problem^ 

The SGS turbulence energy density is defined by the differ- 
ence between the resolved kinetic energy and the filtered kinetic 
energy: 

1 COCO CO 1 ry 1 

K sgs := ~{gu ■ u) A - -g\u\ = --trr sgs , (5) 

where tr T sgs = t„ is the trace of the SGS turbulence stress tensor. 
One can see that the trace of T sgs gives rise to the term -|VA" sgs 
on the right-hand side of Eq. (0. This term can be absorbed into 
the pressure gradient if the thermal pressure P is replaced by the 
effective pressure 

Pea =P+^K sgs =P-itrr sgs . (6) 

The relative contribution of the turbulent pressure P sgs = §A" sgs 
compared to the thermal pressure P is characterized by the SGS 
turbulence Mach number A1 sgs = (IK^/pc 2 ) 112 , where c s is 
the thermal speed of sound. A1 S gs depends on the temperature 
of the fluid and the cuto ff scale A. The dep endence on A is in- 
vestigated in Section l4~2llJoung et al.l (12009) define the turbulent 
pressure by P sgs — (y - 1)K S&S , where y is the adiabatic coeffi- 
cient of the gas. We emphasize that, except for y = 5/3, this 
definition is inconsistent with the decomposition of the fluid- 
dyna mical equations, wh ich fixes the coefficient to be 2/3 (see 
also IChandrasekharl 1951b . This is reasonable because the turbu- 
lent pressure is solely a property of the turbulent flow of a gas at 
a certain length scale, whereas y is a microscopic property of the 
gas that is related to thermal motions of the atoms or molecules. 

The SGS turbulence energy is an intermediate reservoir of 
energy that exchanges energy with the resolved flow and loses 
energy by dissipation into heat. For the computation of A" sgs , a 
PDE has to be solved in addition to the filtered equations for the 
resolved gas dynamics: 

d 

-K sgs + V • (uK sgs ) = T + X - p{e + A) + £. (7) 

While X = TijS ij is the rate of SGS turbulence energy production 
by the turbulent cascade through the cutoff scale A (also called 
the turbulence energy flux) and pe is the viscous dissipation rate 
smoothed over A, effects caused by SGS fluctuations of the grav- 
itational potential and the thermal pressure are given by F and 
pA, respectively. The term D accounts for SGS transport effects. 

3 Althernatively, {qu ® k)a can be expressed in terms of higher-order 
moments. But this merely shifts the closure problem to the higher-order 
moments. 



We refer to ISchmidt et all (l2006al) . Eqs. (33)-(37), for the ex- 
act definitions of these terms. For our purpose it is sufficient to 
discuss the closures of these terms, which are approximations in 
terms of the numerically resolved variables and A"sg S . 

To compute the SGS turbulence stress tensor we propose 
the following closure for the highly compressible regime: 

Tij = 2C 1 A(2 £ >K sgs ) 1/2 S i * - 2C 2 K sgs p^ 2 - |(1 - C 2 )K sgs S u . 

(8) 

where |V ® u\ := (lui^u^) 1 ^ 2 is the norm of the resolved ve- 
locity derivative, S*- = 5y - \ddij is the trace-free part of 

S ij = k( u i,j + u j,d< an d d = iijj. While the first term in Eq. (0 
corresponds to the eddy-viscosity closure that is commonly used 
in incomp ressible LES, the secon d, non-linear term was inves- 
tigated by Wo odward et al.l ([2006) for transonic decaying turbu- 
lence. The standard eddy-viscosity closure follows if C 2 = 0. In 
general, the linear eddy-viscosity term dominates if (K^/p) 1 ^ 2 
is small compared to A\S*\ < A|V ® u\. On the other hand, for 
strong turbulence intensity, i. e., (K sgs /p) 1 ^ 2 > A|V ® u\, the 
non-linear term contributes significantly. This particularly ap- 
plies to intermittent events in supersonic turbulence, for which 
A|V<8>h| > c s . In moderately compressible turbulence, non-linear 
contributions affect the high-intermittency tails of the turbulent 
energy distribution. Independent of the values of C\ and C 2 , 
t„ = -2/T sgs , as required by the identity ©. We denote the trace- 
free part of the SGS turbulence stress tensor by t*.. The verifica- 
tion of the generalized closure <j8j for Ty and the determination 
of the coefficients C\ and C 2 for supersonic turbulence is the key 
to the computation of the turbulent pressure P sgs = | A^, as A" sgs 
first and foremost depends on the production rate £ = TySy in 
Eq. 0. 

Due to the microscopic viscosity v of the fluid, the viscous 

CO 

stresses cry dissipate kinetic energy on the smallest dynamical 

CO 

length scales ( ~ 77 of the physical flow u. The length scale 77 
is called the Kolmogorov scale. In the filtered momentum equa- 
tion (O, viscous dissipation effects are given by the divergence 

00 

of the filtered tensor <x, 7 = (cr (J ) A . The corresponding rate of 
energy dissipation, filtered on the grid scale, is given by 

ge = {o-ijiiij} A = {2vgS u S u ) A = {vg\S u \ ) A . (9) 

It is important to note that ge + o~ijUij, where cr,y and My are the 
filtered viscous stress tensor and the filtered velocity gradient, 
respectively. 

For fully developed incompressible turbulence, the 
Kolmogorov scale can be related to the Reynolds number: 
rj/L ~ Re 3 ^ 4 , where Re := VL/v for an integral length L and 
characteristic velocity V of the flow. As pointed out at the 
beginning of this section, Re is assumed to be very high in 
astrophysical systems. In this case, rj is much smaller than any 
feasible grid resolution A, and simple scaling arguments show 
that the viscous stre ss term in the filtered momentum equa- 
tion (O is negligible dRopke & Schmidtll2009l) . i. e., \cr\ <K |T sgs |. 
Consequently, the physical energy dissipation occurs entirely 
on subgrid scales ( <k A. As 77 decreases in comparison to 
A, the velocity fluctuations on ever smaller length scales give 
rise to arbitrarily steep velocity gradients, which add up to a 
non-vanishing product of the viscosity times the squared rate of 
strain on the right-hand side of Eq. (]9), regardless of how small 
the viscosity is. This results in a non-zero, asymptotically con- 
stant mean rate of energy dissipation in the limit rj — » (v — > 0), 
which is supported by experimental and numerical evidences 
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(see iFrischl 1 19951: llshihara et al] [2009). We may reasonably 
conjecture that the viscous dissipation tensor is negligible in 
the filtered momentum equation and the energy dissipation rate 
does not vanish in the limit of infinite Reynolds numbers also 
in the case of compressible turbulence. A posteriori tests imply 
that this conjecture is fulfilled for driven supersonic turbulence 
(see Sect. I4.2I >. However, the question of energy dissipation 
in inhomogeneous turbulence is more difficult. For example, 
it is known from boundary layers of terrestrial turbulence that 
viscous effects can affect the flow at relatively large scales near 
a wall. For this reason, the microscopic viscosity cannot be 
neglected in LES of such flows. Although solid walls are not 
encountered in astrophysics, many relevant problems exhibit 
pronounced inhomogeneities, and we cannot entirely exclude 
the possibility that viscous effects might become noticeable on 
resolved length scales in certain cases. 

A closure for e follows from simple dimensional reasoning: 



pe : 



C f 



K: 



3/2 



(10) 



Here, it is assumed that the SGS turbulence energy is dissi- 
pated into heat at a rate proportional to K sgs divided by the time 
scale A(K sgs /py l t 2 . For the pressur e-dilatation te r m pA several 
closu res have been proposed (e. g., ISarkarlfT992t iFurebv et al.l 
119971) . However, applying a priori tests (see Sect.[3]l, we find that 
these closures clearly fail in the case of supersonic turbulence. 
The simplest solution i s to neglect pressure dilatation entirely 
dWoodward et al J 12006). In this article, we also set pA - 0, al- 
though we are aware that pressure-dilatation effects have poten- 
tial significance, particularly, in the case of adiabatic turbulence. 
The transport term in Eq. Q can be mo delled by a gradient- 
diffusion approximation (see Sagaut 2003): 



(ID 



where the SGS turbulent diffusivity is approximated by k s „ s w 
0.65A(p^ sgs ) 1/2 , as shown bv lSchmidt etaD (l2006al) . 

In this work, we assume that self-gravity has no significant 
effects on length scales I < A. This corresponds to the condi- 
tion that the local Jeans length A = c s (7r/Gp)^ 2 , where G is the 
gravitatio nal constant, is sufficiently large compared to the grid 
scale A dTruelove et al.lll99llFederrafh et alJl2010al) . Thus, set- 
ting F = 0, the filtered equations resulting from th e compressible 
Navie r-Stokes equations in the limit of rj <sc A (ISchmidt et al.l 
l2006al) read 

-e + v ■(«£>) =0, (12) 

at 

(qu) +V-{gu®u) = p(g + /) - V(P + P sgs ) + V • t* 



8t 



-E + V • (uE) = -£ + pu ■ (g + f) - V ■ [u(P + P sgs )] 



+ V-(h-t* )-2+pe, 



(14) 



where E = \pu + E\ nt is the sum of the resolved kinetic and 
internal energy density, and -X. accounts for sources and sinks 
of the internal energy due to heating and cooling, respectively. 
Since the resolved fluid dynamics on length scales I > A is un- 
affected by the viscosity of the fluid, the above set of equations 
defines the compressible Euler equations for computational fluid 



dynamics in a physically meaningful and consistent way. These 
equations are supplemented by an equation of state, the SGS 
turbulence energy equation (|7), and the Poisson equation for 
the gravitational potential. The pure compressible Euler equa- 
tions without SGS terms, on the other hand, do not follow from 
the compressible Navier-Stokes equation in the limit of infinite 
Reynolds number. In this case, there is no viscous dissipation at 
all, and, by definition, e vanishes identically. This is a mathemat- 
ical idealization that does not describe turbulent flows in nature. 

As a special case, implicit large eddy simulations (ILES) fol- 
low from the above approach. ILES is the most commonly used 
method in astrophysical fluid dynamics. It is based on two as- 
sumptions, which are usually not stated in the literature. Firstly, 
the discretization of the compressible Euler equations introduces 
a dissipative leading error term T) num in the momentum equa- 
tion (fT2l . Implicitly, this term is assumed to be equivalent to 
the SGS turbulence stress term V ■ r sgs . The second assumption 
in ILES is that u ■ T> num = -pe, i. e., kinetic energy on the re- 
solved scales is directly dissipated into heat at a rate that approx- 
imates the viscous dissipation on unresolved length scales. This 
is referred to as numerical viscosity or numerical dissipation. 
Effectively, the following equations are solved in ILES: 



d_ 

dt 



—q + V ■ (uq) = 0, 
at 

(qu) + V -(qu®u) = p(g+f) - VP + D n 
d 



(15) 



(16) 



-E + V-(uE) = -£+pu-(g + f)-V-(uP). (17) 
at 

Despite the lack of a mathematical justification, ILES serves 
as an approximation to turbulent compressible fluid dynamics 
that has proven its utility in numerous astrophysical applications. 
iBenzi et all ((2008) demonstrate that the ILES approach closely 
reproduces two-point statistics of weakly compressible turbu- 
lence in the inertial subrange in comparison to direct numerical 
simulations that solve the Navier-Stokes equations. In this arti- 
cle, we make use of ILES to compute high-resolution data for 
the explicit verification of SGS closures. In contrast, LES treat 
the energy dissipation explicitly. However, it cannot be avoided 
that numerical schemes for compressible fluid dynamics such 
as the piecewise parabolic method (PPM, IColella & Woodward! 
1 19841) introduce some numerical dissipation. Thus, running an 
LES with an explicit SGS model, there is inevitably a numerical 
dissipation channel that competes with the transfer of energy to 
the subgrid scales and the subsequent dissipation of SGS turbu- 
lence energy into heat. Notwithstanding this caveat, we demon- 
strate in this article that physically sensible predictions can by 
made by an explicit SGS model, which are not possible on the 
basis of ILES. 



(13) 3. Closure verification 



To test different closures for the turbulence energy flux S, we 
apply the method described in Schmidt ~et al.l (1200 6a). The basic 
idea is to use data from ILES of non-selfgravitating isothermal 
and adiabatic turbulence with high numerical resolution and to 
apply an explicit filter to these data on a length scale that is in 
between the forcing and the dissipative range. The applied fil- 
ters are Gaussian with filter lengths Ac = 24A and 32A for 768 3 
and 1024 3 grids, respectively. Although the bottleneck effect has 
some influence on the c hosen length scales ( Schmidt et al.1 2009; 
iFederrath et al.l 1201 Obi) , we show that the results change only 
slightly if the filter length decreases or increases by a factor of 
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two. Moreover, comparing to box filters, the results turn out to 
be rather insensitive to the filter type. 



Table 1. Closure and correlation coefficients for the closures 
shown in Figs.Q~|and[2] 



3.1. Single-coefficient closures for supersonic isothermal 
turbulence 

The turbulence energy flux on the filter scale Aq can be com- 
puted explicitly from the unfiltered numerical data by the for- 
mula 



£a g = -guiu j + 



QUiQUi \_ 



(18) 



where the first factor on the right hand side is the turbulence 
stress tensor on the filter scale (defined analogous to Eq. |4| and 
the second factor is the derivative of the filtered velocity. As a 
shorthand notation, we denote the explicitly filtered quantities 
by an overline, for instance, p = (q)a g , and u, = (£>m,)a g /P- 

The above expression for £a g can be compared to closures. 
For the eddy-viscosity closure, Xa g is given by 



I.f o s) = C 1 A G (2p^A G ) 



1/2 



|5 \ 2 --K Ac d, 



where 



1 /— QU 

2 \ Q 



(19) 



(20) 



is the turbulence energy on length scales smaller than Aq. 
Strictly speaking, K Aa is the turbulence energy for the length 
scales ranging from the grid resolution A to the smoothing length 
Aq of the Gaussian filter. If Ac is suffic iently large compared to 
A, this distinction can be neglected (see Sch midt et al.l l2006a). 
Defining 

Ci/ (c.s) = 2 (ds) + 2^- (21) 



the squared error function of the closure can be written as 



err 2 (d) 



■X 



S Ac + ^ Acfl f-C 1 / <cls) 



d 3 x, 



(22) 



where £a g and K Ag are given by Eqs. dTSl and d20l ). respectively, 
and the volume integral extends over the whole domain < V. The 
minimum of err 2 (Ci) yields the least squares error solution for 
the closure coefficient, 



Ci = 



iy\fW\ 2 <P X 



(23) 



where / (ds) = Ag(2p/Ta g ) 1 ^ 2 |5' | 2 for the eddy-viscosity closure. 
For statistically stationary and isotropic turbulence, the closure 
coefficient C\ is independent of the filter length scale because of 
the local equilibrium of the transfer of turbulence energy in the 
inertial subrange. Thus, the value of C\ inferred form Eq. d23l is 
an approximation to the coefficient of the SGS closure for 2 in 
LES. 

To calculate C\, we use data from two 1024 3 simula- 
tions of supersonic isothermal turbule nce with a root-mean - 
square (RMS) Mach number around 5.5 dFederrath et al.l20 10b). 
Statistically stationary and isotropic turbulence is produced 
by stochastic forcing. Solenoidal (divergence-free) forcing is 
applied in one simulation, while the forcing is compressive 
(rotation-free) in the other simulation. We choose Aq = 32A for 
the filtering of the simulation data. Figs.Q]and[2]show the corre- 
lation between and £a g by means of two-dimensional prob- 
ability density functions. For the eddy-viscosity closure ( fT9l ). 



closure 


c, 


corr[Z 3 2A,£^] 


solenoidal forcing, f = 


1.0, Aims 


* 5.3 


eddy viscosity 


0.102 


0.950 


eddy viscosity (shocks excluded) 


0.055 


0.931 


determinant 


0.803 


0.950 


non-linear 


0.849 


0.991 


compressive forcing, ( = 


0.0, , Mrms -5.6 


eddy viscosity 


0.092 


0.930 


eddy viscosity (shocks excluded) 


0.059 


0.914 


determinant 


0.834 


0.947 


non-linear 


0.833 


0.991 



the correlation is quite good (the spacing of the contour lines 
in Figs. Q]and[2]is logarithmic in the two-dimensional probabil- 
ity density), but there is a problem with negative flux values. The 
values of the closure coefficient C\ following from Eq. ( l23T l are 
listed in Table Q] Also listed are the correlation coefficients 



^Pa g -<Sa g >]K s) -<^ s) >] d 3 * 



(24) 



std[E Ac ]std[E^ s) ] 



where std[ ] denotes the standard deviation and the angle brack- 
ets indicate an average over the whole domain. 

An important question for LES of supersonic turbulence is 
whether shocks can be accommodated in closures for the tur- 
bulence energy. Since the e ddy- viscosity closu re originates from 
incompressible turbulence. iMaier et al.l (l2009l) suggested to set 2 
equal to zero in the vicinity of shock fronts. This should suppress 
the spurious production of SGS turbulence energy by the large 
strain at shock fronts. Thus, we tested whether excluding shocks 
in the computation of the eddy-viscosity closure for the turbu- 
lence energy flux would improve the correlations. However, pan- 
els (b) in Figs.[T]and|2]make clear that such a cutoff deteriorates 
the correlations and implies a significant underestimate of large 
positive fluxes. Although the applied shock detection criterion 
d < -c s /Ag is rather crude, we interpret this trend as an in- 
dication that shocks must not be separated from the supersonic 
turbulent cascade. 

In addition to the conventional eddy-viscosity closure, we 
investigate a clo sure that is based on th e determinante of the ve- 
locity gradient (IWoodward et al.l l2~00li). In this case, the trace- 
free part of the SGS turbulence stress tensor is still given by the 
expression t*j = 2pv sgs S*.. The eddy viscosity, however, does 

not depend on K sgs . It is defined by v sgs = -Ci A 2 |5""| ~ 2 det S*. 
Hence, the turbulence energy flux on the filter scale is given by 



r (cls) 
-A G 



(25) 



for this closure. Woo dward et al.l (UpOl) employ the same 
method to test the correlation between their closure and the tur- 
bulence energy flux as we do. A particularly interesting fea- 
ture of the determinant is that it switches signs and thereby de- 
scribes two different flow topologies. In one case, the determi- 
nant is negative. This corresponds to the forward turbulent cas- 
cade transporting energy from large eddies to smaller eddies. In 
the other case, the flow is contracting in one dimension and ex- 
panding in the other two. Then the determinant is positive, cor- 
responding to a backscattering of energy from small eddies to 
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2000 4000 6000 2000 4000 6000 

**3 2 A £3 2 A 

(c) determinant closure (d) non-linear closure 



Fig. 1. Correlation diagrams for the SGS turbulence energy flux in the case of isothermal supersonic turbulence with solenoidal 
forcing (Al rms ~ 5.3). The applied filter length is 32A. The blue dots indicate the average prediction of the closure for a given value 

Of S 32 A- 



larger eddies. This phenomenon can be explained by the align- 
ment of vortices along a single stretching direction (the "tor- 
nado" topology). While an energy flux of the form ( fl9l > fails to 
describe the reverse cascade, we see in panels (c) of Figs.[T]and|2] 
that the determinant closure yields a good correlation for nega- 
tive energy flux. However, the overall correlation does not sig- 
nificantly improve (see Table [TJ, because of the relatively large 
sc atte r in the forward cascade. 

In IWoodward etafl (|2006), a non-linear expression for the 
turbulence stress tensor is investigated, which depends on the 
full Jacobian V ® u of the velocity: 

Illicit j£ 2 

Tjj - -2CiK sgs — — - - -(1 - Ci)K sgs 6ij. (26) 
|V ®m| 2 3 

Since |V ® u\ - (2u.i£Ui t k) the above expression fulfills the 
identity t„ = -2/T sgs . The corresponding turbulence energy flux 



on the filter length scale Ac is given by 

C = ~ 4CiK ^^hf - 3 (1 - c ^ d - (27) 

Figs.[T|(d) and|2](d) show that the correlation is excellent for the 
above closure, with correlation coefficients above 0.99, as listed 
in Table Q] Like the determinant closure discussed above, the 
trace-free part of the non-linear closure for the SGS turbulence 
stress switches signs and, thus, allows for a backward energy 
cascade . 

3.2. Mixed closure for supersonic isothermal turbulence 

Even though the correlation of the turbulence energy flux is very 
good, the purely non-linear closure d26i i is generally not ade- 
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2000 4000 6000 

(a) eddy- viscosity closure 




2000 4000 6000 

(c) determinant closure 



Fig. 2. Correlation diagrams for isothermal supersonic turbulence 

quate as a model for the turbulence stress tensor for the follow- 
ing reasons. Most importantly, rotation invariance is violated be- 
cause of the antisymmetric part of V ® u. As a consequence, 
spurious turbulence energy would be produced for a uniformly 
rotating fluid. Apart from that, the application of this closure in 
LES of forced turbulence show that the growth of turbulence en- 
ergy during the transition from laminar to turbulent flow is insuf- 
ficient for this closure. Thi s is because of the linear dependence 
on /T sgs Q As pointed out bv lWoodward et all d2006h . a seed term 
has to be included in order to trigger the production of turbu- 
lence energy. If the seed term is constructed from the symmetric 
part of the velocity gradient, then turbulence energy production 
vanishes for a uniformly rotating fluid, and, consequently, the 

4 The eddy- viscosity closure depends on K l s ^ . Writing = |pgL s , 
a factor q ms can be cancelled from the SGS turbulence energy Eq. l[7}. 
This results in an equation for q w with a non-vanishing production rate 
of q sgs even starting from the initial condition q sgs = 0. For the non- 
linear closure, on the other hand, g sgs = is a fixed point of the equation. 




2000 4000 6000 

(b) shocks excluded 




2000 4000 6000 

(d) non-linear closure 



with compressive forcing (yVt lms ~ 5.6) as in Fig.Q] 

problem of rotation invariance is also resolved. IWoodward et al.l 
d2006h consider a linear combination of the closure ( f26b with 
the determinant closure. Because of the relatively large scatter 
of the determinant closure, however, we propose a combination 
of the non-linear closure with the linear eddy-viscosity closure. 
Conceptually, this combination has the advantage that the eddy- 
viscosity closure, which is well established for LES of incom- 
pressible turbulence, follows as a limiting case. 

The least-squared-errors approach can be generalized to a 
mixed closure with two coefficient, C\ and C2. For 



Ci/ (cls) + C 2 g {ck> = ££ s) + ^a g d, (28) 
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2000 4000 6000 2000 4000 6000 

(a) solenoidal, Alrms ~ 5.3, Aq = 32A (b) compressive, Mrms ~ 5.6, Aq = 32A 



Fig. 3. Correlation diagrams for the mixed closures for isothermal supersonic turbulence with different forcing. 



Table 2. Closure and correlation coefficients for the linear com- 
bination of the eddy-viscosity and the non-linear closure. 



M 



■A G 



c, 



Co 



corr[£ Ac ,£^'] 



solenoidal forcing (£ = 1.0), At,ms ~ 5.3 



16A 0.96 
32A 1.48 
64A 2.16 



0.0204 
0.0229 
0.0242 



0.749 
0.723 
0.696 



0.991 
0.991 
0.986 



compressive forcing (f = 0.0), Alrms ~ 5.6 



32A 1.29 0.0189 0.698 



0.991 



the closure coefficients are given by the linear system of equa- 
tions 

(jT d 3 *) Cl + (JV^V^d 3 *) C 2 (29) 



= f / (cls) 




§W 


d 3 jc, 


(30) 




ls) d 3 *) 




jV> s >| 2 d 3 *) C 2 


(31) 


= f 8 icls) 






d 3 x, 


(32) 



where 



/ (cls) = A G (2p^ Ac ) 1 / 2 |5*| 2 , 
Ui,kUj,kSjj 



|V®h| 



(33) 
(34) 



The solutions for C\ and C 2 that are obtained from our numer- 
ical data are listed in Table [2] As one can see, the correlation 
coefficients are about as high as for the purely non-linear closure 
and there is only little variation with the forcing and the filtering 
length scale. For Ac = 64A the ratio of Ac to the integral scale 
L is 8, which is quite small. As a consequence, there might be a 
marginal influence of the forcing. The filter length Aq - 16 A, on 
the other hand is significantly affected by numerical dissipation. 



Table 3. Correlation coefficients for the linear combination of 
the eddy-viscosity and the non-linear closure with C\ = 0.02 
and C 2 = 0.7 for isothermal turbulence and adiabatic turbulence 
at various instants with different Mach numbers. 



t/T 


Mm, <A1 2 4A> 


corr[S 2 4 A ,S!,^] 




isothermal (y = 


1.01) 


9.1 


2.2 0.66 


0.990 




adiabatic (y = 


1.4) 


2.0 


1.3 0.48 


0.981 


3.9 


0.9 0.34 


0.986 


8.0 


0.6 0.26 


0.990 


15.9 


0.5 0.21 


0.990 



The correlation diagrams for the mixed closure are plotted in 
Fig. [3] Although the relation between Erf? and Y32A is slightly 
tilted for negative fluxes, the results are comparable to Q] (d) 
and|2](d). Therefore, we base our SGS model on the mixed non- 
linear closure ([8} with the averaged coefficients C\ = 0.02 and 
C 2 = 0.7. 



3.3. Supplementary tests for different Mach numbers 

For the data listed in Table|2] the average Mach numbers associ- 
ated with the filter scale, Ma g = ((2^a g /pCj) 1 ^ 2 ), assume values 
around the speed of sound. Thus, the question arises whether the 
closure coefficients calculated above are applicable to subsonic 
velocity fluctuations. To test the mixed closure for a different 
Mach number, we calculated the turbulence energy flux S^ s) 
with fixed values C\ = 0.02 and C 2 = 0.7 for dat a from a simu- 
lation of isothermal turbulence with Alrms ~ 2.2 (Schmid t et al.l 
2009). The resulting correlation diagram is plotted in Fig.|4](a). 
There is a small bias to overestimate the turbulence energy flux, 
but the prediction of the mixed closure is still very good. Indeed, 
a correlation coefficient coit[Sg, £q 1s) ] = 0.990 is obtained in 
this case (see TableO. In addition , we investigated da ta from an 
adiabatic turbulence simulation (Schmi dt et al. 1 12007b . in which 
the RMS Mach number gradually decreases with time because 
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of the dissipative heating of the gas. The results are summarized 
in Table [3] and the correlation diagram for the final snapshot of 
the simulation is shown in Fig. 0(b). Our results suggest that the 
closure coefficients are not very sensitive to the Mach number. 
Nevertheless, we cannot exclude that the optimal values of C\ 
and C2 differ significantly if the SGS turbulence Mach number 
Msgs (see Section[2]i is only a tiny fraction of the speed of sound. 
Answering this question is left for future studies. 

3.4. Energy dissipation 

For the turbulence energy on the length scale Ac, which is de- 
fined by Eq. (l20l i. a dynamical equation analogous to Eq. (0 can 
be formulated. Averaging this equation over the whole periodic 
domain and assuming statistical equilibrium, i. e., d,(K& G ) =s 0, 
the following global balance equation is obtained: 

(Z Ao ) - (Pd -Pd) - C^-^-j ^ 0. (35) 

The first term is the mean turbulence e nergy flux, the secon d 
term is the mean pressure dilatation (see Schmidt et al. 2006a), 
and the third term is the mean dissipation rate expressed in terms 
of K& G . Substituting Eq. (TT~8T > for Sa g , yields the coefficient of tur- 
bulence energy dissipation, C e . From the supersonic isothermal 
turbulence data, we find a value C 6 » 1.5, which is somewhat 
higher yet still comparable t o typical valu es calculated for in- 
compressible turbulence (see Sagaut 2006). 

4. Large eddy simulations of forced supersonic 
turbulence 

To investigate statistical properties of the SGS turbulence energy 
and related quantities, we run LES of forced supersonic isother- 
mal turbulence with the SGS model defined in Sects. [2] and [3] 
For the implementation, we use the code Enzo 1 .5 developed by 
the Laboratory for Computational Astrophysics at the University 



of California in San Diego (http://lca.ucsd.edu). In these simula- 
tions, we apply solenoidal, compressive and mixed force fields 
to produce statistically stationary and homogeneous turbulence 
with different RMS Mach numbers (see ISchmidt et ail 1.2006; 
ISchmidt et alJ l2009: Fed errath et al.ll2010bl) . The forcing acts on 
length scales around the integral length L, where L is one half of 
the box size. The autocorrelation time of the force field is given 
by the time scale T = L/V, where the characteristic velocity V 
specifies the magnitude of the turbulent velocity fluctuations on 
the integral scale. The mixture of solenoidal (divergence-free) 
and compressive (rotation-free) modes of the force field is ad- 
justed by means of a Helmholtz decomposition with weighing 
parameter < £ < 1. Purely solenoidal forcing results for 
( = 1. Setting the adiabatic exponent y = 1.001, the energy 
dissipated per integral time is small compared to the internal en- 
ergy for Atrms up to about 10. For this reason, the gas is pseudo- 
isothermal. This approximate treatment of isothermality enables 
us to monitor energy conservation. With our implementation of 
the SGS model, the sum of resolved kinetic energy, SGS turbu- 
lence energy, and internal energy minus the power of the forcing 
integrated over time is conserved for the whole computational 
domain to a relative precision better than 10~ 8 . The fraction of 
computational time consumed by the SGS model is in the per- 
cent range. 

4.1. Correlations with resolved flow quantities and the 
effective pressure 

As an example, Fig. [5] shows a visualization of K sgs prepared 
from an LES with 512 3 grid cells. The parameters of this simu- 
lation were chosen to match the ILES with solenoidal forcing 
in Section I3.ll The RMS Mach number of the flow is about 
5.5 in the statistically stationary regime. In the reddish regions, 
A" sgs is higher than the spatial mean, while it is lower in the 
bluish regions. For comparison, Fig|6]shows the local denstrophy 

Q 1/2 = I V X [q 1/2 u)\ , which is an indicator of compressible 
turbulent velocity fluctuations (Kritsuk et aD l2007h . It appears 
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Fig. 5. Visualization of the SGS turbulence en- 
ergy density K sgs in a 512 3 LES with solenoidal 
forcing. 



that high SGS turbulence energy is concentrated in regions of 
intense denstrophy. On average, K sgs ~ 0.1A 2 Qi/2 for large den- 
strophy values, as one can see in the correlation diagram of /T sgs 
vs. A 2 Qi/2 in Fig. |7] (a). The same relation is found for com- 
pressive forcing (see panel (b) of Fig. |7j. Nevertheless, the local 
values of K sgs and A 2 Qi/2 deviate substantially from the aver- 
age relation. This is a consequence of the various processes con- 
tributing to the SGS dynamics, which are not fully encompassed 
by the derivative of the resolved velocity field. For this reason, 
derived quantities such as the rate of strain or the denstrophy are 
only of limited utility to estimate effects of turbulence on unre- 
solved length scales. Since f sgs = §^ S gs, this applies also to the 
turbulent pressure. 

The phase diagrams of the effective pressure (0 vs. the mass 
density are plotted in Fig. [8]for both LES. One can see that the 
average of the effective pressure for a given mass density closely 
follows the isothermal relation P oc p. This is because the mean 
turbulent pressure is small compared to the thermal pressure for 
the resolution A = L/256 (see Section l4~2l i. Locally, however, 
the intermittency of turbulent velocity fluctuations can give rise 
to an effective pressure that exceeds the thermal pressure by one 
order of magnitude. For this reason, the contribution of the tur- 
bulent pressure P sgs is locally not negligible. This effect becomes 
stronger as the cutoff scale A increases in comparison to the in- 
tegral scale of turbulence. 

In Sect. [2] we argue that the viscous stress term in the 
filtered momentum equation (0 vanishes in the limit of infi- 
nite Reynolds number and the rate of energy dissipation on the 
grid scale, e, is determined by the SGS turbulence energy (see 
Eq. [Tol l. In contrast, an extrapolation of the expression for the 
microscopic dissipation rate on length s cales I ~ r\ to the grid 
scale A was proposed bv lPan et alJ ( 20091) : 



The grid-scale viscosity va = const, in the above expression 
is treated as a constant coefficient that is determined by the 
mean numerical dissipation of PPM. Defining the compressible 
Reynolds number of the resolved flow by Re a = 2L 2 (\S *| 2 )/M 2 ms , 
where u llns is the root mean square velocity the viscosity can be 
evaluated from va = VL/Re&. The problem with this approach 
is that the viscosity on the grid scale, which corresponds to the 
SGS eddy-viscosity, cannot be assumed to be constant. 

Neglecting diffusion, compressibility and the non-linear 
term in the SGS turbulence stress ©, the equilibrium between 
production and dissipation of SGS turbulence energy in Eq. ^ 
implies K sgs ~ (Ci/C e )£>A 2 |S*| 2 . Hence, e ~ (A/C e ) 2 (d|S*|) 3 
according to Eq. (fTOl . We emphasize that a relation of the form 
e ~ A 2 !^*! 3 follows from any co mmon SGS m odel under the 
assumption of local equilibrium (Sagaul l2006l) . Comparing to 
Eq. (|36T >. we see that that va ~ A 2 |5*|, which is not a constant. 
This is a consequence of the fact that ge + o~ijiiij 00 \S*\ 2 , as ex- 
plained in Sect [2] The discrepancy becomes apparent in Fig. [9] 
which shows the correlation diagrams of the rate of energy dis- 
sipation calculated via Eq. d36l i vs. e following from the SGS 
model. Toward low values of e, we find an average relation close 
to |5*| 2 oc e 2 ^ 3 , which is just the relation that follows from the 
above estimate of the equilibrium dissipation rate. This behavior 
is reasonable because the contribution of the non-linear term in 
the closure (0, which is neglected in the estimate, is relatively 
small for low values of K sgs (corresponding to low energy dis- 
sipation). Moreover, the unresolved velocity fluctuations tend to 
be small compared to the speed of sound in this limit, which cor- 
responds to low compressibility. Consequently, the results from 
the LES support the theoretical arguments against Eq. (l36l i as an 
approximation to the dissipation rate on the grid scale. Although 
Eqs. ( [Tol l and (l36l i yield about the same mean dissipation rate, 



ge = gv A \S' 



(36) 



5 See lSchmidt et alJ d2009ri . Here, we replace u) 2 ^ by <|5 *| 2 > = {or + 
d 2 ) for consistency with Eq. i36\ . 
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Fig. 6. Visualization of the denstrophy ili/2 for 
the same snapshot as in Fig. [5] 
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log 10 (A 2 Xl 1/2 /puV 2 ) 
(a) solenoidal 
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log |0 (A 2 ni /2 /p ^ 2 ) 
(b) compressive 



Fig. 7. Correlation diagrams of the SGS turbulence energy vs. the denstrophy, normalized by the cutoff scale A, for 512 3 LES with 
solenoidal and compressive forcing. The contours are logarithmic. The average relation between both quantities is indicated by the 
dotted lines, and the dashed line shows the relation /T sgs ~ 0.1A 2 Qi/2. 



the former determines the local rate of energy dissipation on the 
footing of a physically well motivated scale-separation of fluid 
dynamics, while the latter is based on a putative analogy between 
the numerical and the microscopic viscosity. 



4.2. Dependence on the cutoff scale 

The scaling of the turbulent velocity fluctuations in su- 
personic hydrodynamic turbulence has been inferred from 
energy spectrum f unctions and structure functions (e. g., 
Krits uk et al J 120071: ISchmidt et all 120081 120091: iFederrath et ail 
201 Obi iPrice & Federrathll2010h . The pure velocity scaling in 
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log J0 (p/p ) log J0 (p/p ) 
(a) solenoidal (b) compressive 

Fig. 8. Phase diagrams of the effective pressure defined by Eq. (O vs. the mass density for 5 12 3 LES with solenoidal and compressive 
forcing. The contours are logarithmic. The averages of the SGS turbulence energy for particular values of the denstrophy are 
indicated by the dotted lines. 
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Fig. 9. Correlation diagrams of the normalized rate of energy dissipation defined by Eq. ( l36l l. where va assumes a constant value 
that is given by the numerical Reynolds number, vs. the rate of energy dissipation ( fTUb that is predicted by the SGS model. The 
averages of expression (l36l l for given values of pe sgs are indicated by the dotted lines. 



the supersonic regime is stiffer than Kolmogorov scaling, and 
it appears that the scaling expone nt depends on the forcing. For 
example. iFederrath et al.1 (l2010bl) find indices of the turbulence 
energy spectra f3 = -1 .86 + 0.05 and -1 .94 + 0.05 for solenoidal 
and compressive forcing, respectively. For incompressible tur- 
bulence, p = -5/3. Velocity variables with fractional mass- 
weighing, in particular p l ^u, exhibit similar scaling laws, which 



can be interpreted as an indication of universality (Kritsu ket al.l 
120071 ISchrmdUtal]l2008l) . 

The SGS turbulence energy is given by the fluctuations of 
the velocity and density fields on length scales ( < A, as defined 
by Eqs. (|4]l and (0. However, there is no obvious relation to the 
known scaling laws for turbulence, because the decomposition 
of the fluid dynamical variables cannot be related to the two- 
point statistics (structure functions) or the Fourier modes (energy 




spectra) in a straight-forward manner. To determine the scaling 
of SGS turbulence as a function of A, we run several LES with A 
ranging from L/256 to L/32. The mean values of the RMS Mach 
number and the SGS turbulence Mach number are plotted as 
functions of time in Fig. [10] The flow approache s a statistically 
stationary state after about 2 integral time scales dSchmidt et alj 
2009; Fed errath et ai.ll20" l0b), for which Ai rms settles at values 
between 5 and 6. The temporal variation of Ai rms is caused by 
the stochastic forcing. As expected, (Ml os ) 1/2 decreases with the 
cutoff scale. Averaging the spatial means from t — 2T to 107\ we 
find the time-averaged mean values listed in Table[4] As one can 
see in Fig.QT](a), the ti me averages of (M 2 g J 1 ^ 2 closely follow 
power laws, 

(Msgs)'' 2 ^ < 37 > 

with a M = 0.475 + 0.004 for solenoidal and 0.451 + 0.026 for 
compressive forcing. 

The behaviour of the mean SGS turbulence energy is sim- 
ilar, although the intermittent fluctuations of (K sgs ) are more 
pronounced in comparison to the mean SGS turbulence Mach 
number (see top panels of Fig. fT2l . The higher degree of inter- 
mittency stems from the mass density that is included in K sgs . 
For compressive forcing, (K sgs ) is systematically lower in com- 
parison to the LES with solenoidal forcing. This indicates that 
the total amount of energy in the turbulent structures on a given 
length scale is smaller in the compressive forcing case. The ra- 
tio of the mean values of K sgs for compressive and solenoidal 
forcing in Table [4] approximately agree with the ratio 0.38 that 
is inferred from the filtered high-resolution data. On the other 
hand, the scaling laws 

(U«A«', (38) 



are nearly the same for solenodial and compressive forcing (see 
Fig. E] (b)). We find the slopes a K = 0.799 ± 0.009 and 
0.769 + 0.029, which agree within the error bars. This result is 
remarkable, because it suggests that the scaling properties of tur- 
bulence on small length scales are independent of the forcing. 

As can be seen in the bottom panels of Fig. [12] the above 
scaling law of (K sgs ) results in a mean dissipation rate (pe) that 
is independent of the cutoff scale, which is an essential property 
of the energy dissipation predicted by the SGS model. The time- 
averaged mean values are listed in Table [4] The significantly 
lower mean dissipation rate in the case of compressive forc- 
in g is consistent with the energy spectra of p l ^ 3 u (see Fig. A.l 
in lFederrath et al. 2010b). This mass-weigh ted velocity variabl e 
is related to the energy dissipation rate dKritsuk et al.l 12007b . 
Moreover, Fig. [13] (right panel) shows that the growth of the 
mean internal energy in time becomes smaller as the weighing 
parameter £ decreases from 1 (solenodial forcing) to (compres- 
sive forcing). After subtracting the c ontribution from num eri- 
cally resolved compression effects (see lSchmidt et al.ll2006l) . we 
find that, independent of the cutoff length A, about 3/4 of the 
change of the internal energy stems from SGS turbulence en- 
ergy dissipation. The remainder is caused by numerical dissipa- 
tion. This does not imply that the total rate of energy dissipation 
is much higher in LES compared to ILES, because the total en- 
ergy dissipation is always determined by the energy injection 
due to the forcing. In conclusion, the greater part of kinetic en- 
ergy is dissipated through the SGS turbulence energy reservoir 
at a scale-free rate. 

To quantify the relative importance of large values of Ai sgs , 
we determine the volume fractions of cells with an SGS turbu- 
lence Mach number greater than a particular value. This fraction 
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Fig. 11. Scaling laws for the mean SGS turbulence Mach number (a) and energy (b) as functions of the numerical resolution A. 

Table 4. Time-averaged spatial mean values of various quantities and their standard deviations from the averages for different 
numerical resolutions. 



N 


A/L 


Aims 






solenodial forcing = 1) 


64 


1/32 


5.38 


1.107 +0.053 0.0726 ± 0.0055 


1.236 ±0.141 


128 


1/64 


5.50 


0.787 ± 0.030 0.0407 ± 0.0025 


1.230 ±0.113 


256 


1/128 


5.55 


0.578 ± 0.022 0.0236 ± 0.0013 


1.213 ±0.098 


512 


1/256 


5.52 


0.412 ±0.012 0.0138 ±0.0012 


1.219 ±0.159 


compressive forcing = 0) 


64 


1/32 


5.29 


1.353 ± 0.049 0.0235 ± 0.0044 


0.253 ± 0.072 


128 


1/64 


5.43 


1.040 ±0.041 0.0148 ±0.0018 


0.286 ± 0.047 


256 


1/128 


5.57 


0.767 ± 0.029 0.0086 ± 0.0013 


0.293 ±0.061 


512 


1/256 


5.86 


0.528 ± 0.023 0.0048 ± 0.0008 


0.292 ± 0.067 



is given by 1 - cdf(At sgs ), where cdf(A1 sgs ) is the cumulative 
distribution function of M sgs - In Figure [14] the resulting func- 
tions are plotted for the LES with different cutoff lengths. As 
expected, the fraction with Af sgs > 1 decreases with the cutoff 
length A. However, the tails toward high AL, gs demonstrate that 
even at relatively high resolution there are supersonic velocity 
fluctuations on unresolved length scales, and the corresponding 
turbulent pressure decreases only little with the cutoff scale. 

For the lowest-resolution LES, we can compare the dis- 
tribution of Alsgs to the distribution inferred from the corre- 
sponding filtered 1024 3 data (see Section |3~TT >. The filter length 
Ac = 16A = L/32 is equivalent to the cutoff length in the 64 3 
LES. Choosing yet a lower resolution of the LES, correspond- 
ing to a larger filter length for the ILES, turned out not to be 
feasible. Even for the LES with A/L = 1 /32, the forcing range 
and the range of length scales that are directly affected by nu- 
merical dissipation overlap. For the filtering of the ILES, on the 
other hand, the filter length cannot be lowered (corresponding to 
a higher resolution of the LES), because the dynamical range of 
fluctuations between the grid scale and the filter length would 
become insufficient and the numerical smoothing would be too 
strong. Nevertheless, Fig. [14] (a) demonstrates that the distribu- 
tions agree remarkably well is the case of solenoidal forcing. For 
compressive forcing, there are larger discrepancies. However, 
given that Gaussian filtering corresponds only roughly to the im- 
plicit filter in an LES and that the SGS model is based on various 
approximations, the match is quite satisfactory. The larger devi- 
ations in the case of compressive forcing suggest that it is not 
possible to calibrate the SGS model coefficients in such a way 
that an optimal match is obtained both for solenoidal and for 
compressive forcing at the same time. The different shape of the 



distribution that is obtained from the high-resolution simulation 
with compressive forcing points toward a missing physical ef- 
fect such as the pressure-dilatation, which is entirely neglected in 
our SGS model. Anyhow, purely compressive forcing is a limit- 
ing case. In nature, some mixture of solenoidal and compressive 
forcing is more likely to occur. In Fig. [13] (right panel), we com- 
pare the distributions of A1 sgs for force fields with £ varying from 
1 (solenoidal) to (compressive). High SGS turbulence Mach 
numbers become more frequent as the contribution of compres- 
sive forcing modes increases. 

4.3. Dependence on the Mach number 

At fixed resolution, the SGS turbulence energy increases with the 
resolved kinetic energy of the flow. For isothermal turbulence, 
this also implies an increase of the SGS turbulence Mach number 
with rising RMS Mach number. To investigate this dependence, 
we varied the magnitude of mixed forcing with £ = 2/3 and 
f = 1 /3. In the case ( - 1 /2, a forcing field with two solenoidal 
and one longitudinal de grees of freedom is obtained. Fig. 8 in 
iFederrafh et al.l (1201 0b) demonstrates that the ratio of the en- 
ergy that is contained in transversal and longitudi nal modes ap- 
proac hes a constant value for £ > 1/2 (see also Rritsu ket al.l 
2010). Correspondingly, the cumulative distributions of AL, gs 
plotted in Fig. [13] (b) show that there is almost no difference 
between forcing with f = 2/3 and purely solenoidal forcing 
= 1). On the other hand, one can see that there is a noticeable 
influence of compressive modes for £ = 1/3, but the distribution 
differs from the purely compressive case (£ = 0). 

The dependence of the RMS SGS turbulence Mach number 
on the forcing amplitude is shown in the top panels of Fig.fTBIfor 
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Fig. 12. Temporal evolution of the SGS turbulence energy (top) and the dissipation rate (bottom) for solenoidal (left column) and 
compressive forcing (right column). The cutoff length A decreases from L/32 (light colour) to L/256 (full colour). 



256 LES. The corresponding characteristic Mach numbers and 
time-averaged statistics are listed in Table [5] Independent of the 
Mach number, the ratio of (Ml gs ) 1 ^ 2 to A1 lms is nearly equal for 
£ = 2/3 and 1/3. This ratio is 0.102 for purely solenoidal forc- 
ing = 1), and 0.138 for purely compressive forcing (£ = 0). 
Consequently, <A^ 2 os )^ 2 /A1 rms is mostly determined by the grid 
scale A, except for small As one can see in the middle pan- 
els of Fig. [15] the normalized mean SGS turbulence energy, 
(Ksgs) / puV 2 , is about the same for the different Mach numbers, 
with a weak trend to decrease toward low Mach numbers (see 
also Table |5)- The same behavior is found for the mean dissi- 
pation rate (Fig. [15] bottom panels), which further supports the 
validity of the SGS model in the supersonic regime. Following 
the trend discussed in Sect 14.21 there is clearly an influence of 
the mixture of solenoidal and compressive modes in the forcing. 

Fig. [TBI shows the distributions of A1 sgs for £ = 2/3 (a) and 
1/3 (b), as explained in Section l4~2l For both forcing types, the 
volume fractions with At sgs > 1 increase with the forcing mag- 
nitude. For RMS Mach numbers greater than 5, supersonic tur- 
bulent velocity fluctuations at the cutoff scale fill more then 10 % 
of the total volume. If the ratio between the integral scale of tur- 
bulence and the cutoff scale is smaller, the volume filling factor 
increases further. 



5. Conclusion 

Formulating a mixed closure for the flux of energy from the nu- 
merically resolved to the unresolved scales, we have generalized 
the subgrid-scale (SGS) turbulence energy model to the regime 
of highly compressible tu rbulence. This closure is based on ideas 
of Wood ward et alJ (120061) and features a non-linear term in ad- 
dition to a linear eddy-viscosity term. In general, the turbulence 
energy cascade is an important source of SGS turbulence energy 



production in turbulent flows, and it should not be neglected 
even if other sub-resolution sources a re emphasized in partic- 
ular astrophysical applic ations (e. g., Scannapieco & Briiggen 
2008; Joung et al. 2009). Our proposed closure for the trans- 
fer of energy by the turbulent cascade complements theses 
models. We verified this closure by means of explicit filtering 
of high-resolution data from various simulations of supersonic 
isotherm al and adi abatic turbulence (Schmidt et al. 2007, 2009; 
Federrath et a U2010bl) . Tests in large eddy simulations of forced 
supersonic turbulence show that the SGS model meets several 
important requirements that should be satisfied by any sound 
SGS model: 

- For statistically stationary turbulence, an equilibrium be- 
tween the production and the dissipation of SGS turbulence 
is reached. The mean SGS turbulence energy depends on the 
grid scale via a power law. 

- The SGS turbulence Mach number, which specifies the im- 
portance of the turbulent pressure on the grid scale relative 
to the thermal pressure of the gas, depends linearly on the 
RMS Mach number of the resolved turbulence. 

- The SGS turbulence energy dissipation is independent of the 
grid scale. 

Forced turbulence simulations in a periodic box are most suit- 
able to test the properties listed above, because of the well de- 
fined statistics of the isotropic, homogeneous and stationary tur- 
bulence that is produced. In addition to these properties, we 
found a dependence of the SGS turbulence energy and the rate 
of energy dissipation on the mixture of solenoidal (divergence- 
free) and compressive (rotation-free) forcing modes. However, 
the scaling laws for the SGS turbulence energy are very similar 
for solenoidal and compressive forcing. This is an indication that 
the SGS model describes the dynamics in the inertial subrange, 
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Fig. 13. Left panel: Time evolution of the mean internal energy for forcing with varying £ and about the same RMS Mach number. 
Right panel: Volume fractions of zones, in which the SGS turbulence Mach number is greater than a certain value for the same 
forcing parameters as in the left panel (the curves for ( = 1 and £ = 2/3 almost coincide). 



l-cdf(Al, s ) l-cdf(M sgs ) 




(a) solenoidal (b) compressive 



Fig. 14. Volume fractions of zones, in which the SGS turbulence Mach number is greater than a certain value for different numerical 
resolutions. The dashed lines follows from the filtering of 1024 3 ILES data with filter length 16A. This corresponds to the 64 3 LES, 
for which A/L = 1 /32. 



Table 5. Time-averaged spatial mean values of various quantities and their standard deviations from the averages for different 
forcing magnitudes (defined by Ma = V/cq) and mixtures of solenoidal and compressive modes. 



Ma 


Aims 






(K m 


)lpoV l 


(.L/V'Xe) 


f = 2/3 


3.33 


2.91 


0.293 +0.010 


0.101 


0.0193 


± 0.0016 


0.933 


±0.109 


6.67 


5.38 


0.559 ± 0.021 


0.104 


0.0223 


±0.0016 


1.125 


±0.123 


10.00 


7.54 


0.778 ± 0.044 


0.103 


0.0229 


±0.0018 


1.150 


±0.142 


16.67 


10.49 


1.097 ±0.117 


0.105 


0.0232 


± 0.0017 


1.168 


±0.130 


f =1/3 


4.25 


2.96 


0.320 ± 0.009 


0.108 


0.0129 


±0.0011 


0.535 


± 0.064 


8.50 


5.73 


0.636 ± 0.020 


0.111 


0.0149 


± 0.0015 


0.641 


± 0.089 


12.75 


8.14 


0.896 ± 0.052 


0.110 


0.0153 


± 0.0017 


0.666 


±0.101 


21.25 


11.81 


1.268 ±0.123 


0.107 


0.0152 


± 0.0015 


0.654 


± 0.094 



although the length scales close to the cutoff scale are affected by 
numerical dissipation. The differences in the mean values result 
from the substantial diffe rences in the turbulent flow s tructure on 
larger scales (see lFederrath et al.l l2008. 2009. l2010bl) . 

The implementation of the SGS model into a fluid-dynamical 
code such as Enzo is rather straightforward. The SGS turbulence 
energy can be treated as a passive scalar with various source 
terms. To evaluate the mixed closure for the turbulence energy 
cascade, derivatives of the resolved velocity field are computed 
by means of centered differences. Care must be taken to en- 
sure energy conservation, particularly, if the net change of the 



SGS turbulence energy in a certain grid cell exhausts the avail- 
able energy over a time step or, vice-versa, if too much energy 
is drained from the resolved scales. However, these are excep- 
tions that can be handled numerically. The effective increase in 
computing time is less than 10 %. Moreover, the hydrodynamic 
coupling of the SGS model to the resolved flow introduces a re- 
duction of the bottleneck effect in the turbulence energy spectra 
dWoodward et al.ll2006l: [Schmidtll20 1 Ot) . 

In our large eddy simulations, we find correlations between 
resolved flow quantities such as the rate of strain or the denstro- 
phy and SGS quantities, but the scatter is large. This is problem- 
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Fig. 15. Temporal evolution of the SGS turbulence Mach number (top) and energy (middle) and the dissipation rate (bottom) for 
f = 2/3 (left column) and ( = 1/3 (right column). The different lines in each plot correspond to different forcing magnitudes (see 
Tabled, which are specified by the values of the characteristic Mach number Ma = V/cq (cq is the initial speed of sound). 



atic if one intends to estimate unresolved flow properties on the 
basis of such correlations. In particular, this applies to the calcu- 
lation of the dissipation rate from a c onstant numerica l viscosity 
and the rate of strain, as proposed bv lPan et all J2009). We have 
shown that the assumption of a constant numerical dissipation 
coefficient is inconsistent with the equilibrium relation between 
the dissipation rate and the rate of strain on the grid scale in the 
limit of large Reynolds numbers. This relation, which follows 
from the SGS turbulence energy model with the linear eddy- 
viscosity closure (and also from the Smagorinsky model), is ver- 
ified by our LES data for low turbulence intensity, while devia- 
tions become apparent for strong turbulent dissipation. This can 
be understood as a consequence of the non-linear term in the 
closure for the turbulence energy flux. Also the estimate of tur- 
bulent pressure effects on the basis of the rate of stra in and the 
vorticity of the resolved flow that is put forward by IZhu et all 
d2010l) is incomplete because they do not distinguish between 



the contribution from the resolved flow and from the subgrid 
scales. The predictions from both approaches with regard to tur- 
bulen ce in the intergalactic medium are compared in ongoing 
work (Iap ichino et al.ll20loh . 

From the probability distributions of the SGS turbulence 
Mach number, it follows that the turbulent pressure locally ex- 
ceeds the thermal pressure even at moderate RMS Mach num- 
bers and for relatively small grid scales. Sin ce the grid scale in 
contemporary galac tic disk simulations (e. g jAgertz et al.l2 009: 
Tasker & Tan 2009) is close to molecular cloud scales (a few pc), 
unresolved supersonic velocity fluctuations are quite likely and 
the turbulent pressure plays an important role. This has implica- 
tions for the treatment of collapsing gas regions. The criterion 
for gravitational stability, which influences the grid resolution 
in adaptive mesh refinement simulations and controls the pro- 
duction of sink particles to capture the collapsing gas, is usu- 
ally based on the thermal Jeans mass (among other criteria; see 
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Fig. 16. Volume fractions of zones, in which the SGS turbulence Mach number is greater than a certain value for different RMS 
Mach numbers and mixtures of solenoidal and compressive modes, as in Fig. IT31 



iFederrath et al.ll2010al) . To account for the effects of turbulence 
below the grid scale, we suggest to include the turbulent pressure 
in the definition of the Jeans mass, in analogy to the magnetic 
pressure in self-gravitating MHD turbulence. 

To model the fragmentation below the grid scale in more de- 
tail, a possible approach is based on the assumption that the lo- 
cal star formation efficiency is regulated by turbulence on the 
grid scale. Then the star formation rate can be parameterized in 
terms of the tur bulent Mach number that is calculated from the 
SGS m odel (see lKrumholz & McKed2005l:lPadoan & Nordluridl 
2009). On the other hand, star formation acts back on the 
SGS turbul e nce e nergy via stellar feedback. As suggested by 
I Joung et all (|2009), a stellar feedback term can be included into 
the SGS turbulence energy equation. Statistically, we have 



d 

dt 6 



£ + X* - pe, 



where X* oc pe+frg accounts for the energy injection per unit 
mass, e*, by supernovae. The associated time scale is the free- 
fall time scale Ts = [3n/(32Gp)] 1 ^ 2 , which is the fundamental 
time scale of star formation. Neglecting the fluctuations of the 
gas density and setting the mean production rate S ~ V 2 /T, 
where V and T are the typical velocity and the turn-over time 
scale, respectively, of the resolved turbulent flow, it follows that 
the turbulent pressure in equilibrium is of the order 



— +C* — 

T r ff 



2/3 



(39) 



Kolmogorov scaling becomes manifest in the factor A 2 ^ 3 in 
Eq. d39l . For highly compressible turbulence, however, the scal- 
ing of the SGS turbulence energy deviates from the Kolmogorov 
law (see Sect. 14. 2\ . Depending on the ratios V 2 /e* and T/ts, the 
production of SGS turbulence energy by the turbulent cascade or 
by supernovae dominates. The model of Joung et al. (200§|) fol- 
lows in the limit (V 2 /e*)(Tff /T) <k 1. In general, shear instabil- 
ities, gravitational instabilities, cooling instabilities, etc. above 
the grid scale feed energy to smaller scales. For disk galaxies, a 
simple estimate can be obtained from the v elocity dispersion of 
atomic hydrogen, which is about 10km s'. lAgertz et alJ (T2009) 
show that gravitational instabilities grow on length scales rang- 
ing from 0.1 to about 2kpc. Setting V = 10 6 cms _1 and assum- 
ing L > 0.1 kpc » 3 x 10 20 cm, the turbulence energy flux to 
smaller length scales is V 2 /T = V 3 /L < 0.003 erg g -1 s _1 . On 
the other hand, <?* « 4 x 10 48 T/ergMg 1 « 2 x 10 15 ergg~' and 



C* » 0.025 imp ly ~ 0.03/7 ("/ 1 cnr 3 ) 1/2 erg g" 1 s _1 (see 
iJoung et ai1l2009l) . Since the efficiency of the energy transfer 
from sup ernova blast waves to the i nterstellar gas is roughly 
77 ^ 0.1 (Ma c Low & Klessenl 120041) . the energy injection by 
stellar feedback is comparable to the turbulence energy flux for 
atomic hydrogen with density n ~ lcirT 3 . The dependence of 
2* on the gas density implies a greater contribution form super- 
nova feedback in the cold gas phase, but the above estimate does 
not account for the intermittency of turbulent velocity fluctua- 
tions, which entails large deviations from th e mean. As a conse- 
quence, the assumption of Joung et al. (2009) to consider the en- 
ergy injection by supernova as the main source of the turbulent 
pressure is marginally fulfilled in cosmological simulations, in 
which the internal structure of galaxies is very poorly resolved. 
In galactic-scale simulations with high resolution, on the other 
hand, turbulence is not uniformly produced, and including the 
turbulence energy cascade improves the description of numeri- 
cally unresolved processes. In particular, it will be useful to at- 
tempt a further generalization of the SGS model to multi-phase 
turbulence. A very s imple ansatz has recently been presented by 
iMurante et al.l (120101) . A complete SGS model that accounts for 
a warm and a cold gas phase is presently under development 
dBraun & Schmidtl2010l) . 

Including stellar feedback and cooling into our SGS 
model will be of further utility for the numerical treat- 
ment of turbulence in the intergalactic medium (see 
Scringel & Hernquist 2003), where turb ulence is produced 
by different processes (see, for instance, ICen & Ostrikerl l 1999; 
Subraman ian et al.ll2006t iRvu et alj|2008t llapichino et al.l l2008, 
2010). lOppenheimer & Pavel (120091) show that a significant 
amount of the line broadening of O VI in cosmological simula- 
tions stems from numerically unresolved turbulence. They apply 
a heuristic model in the postprocessing of the simulation data. 
Using an SGS model, on the other hand, the effect on the line 
broadening can be computed on the fly. Moreover, metals are 
mixed into the intergalactic medium by turbulence that is driven 
by galactic outflows. SGS turbulence enhances the turbulent 
mixing. Following an approach that is quite similar to the treat- 
ment of stellar feedback in galaxy simulations, our SGS model 
can be used in combination with phenomenological models for 
super nova-driven outflows dJoung et al.1 120091 : lEvoli & Ferraral 
2010). In both cases, the use of adaptive mesh refinement is 
mandatory to achieve a sufficient dynamical range. Therefore, 
an essential objective for future work will be to incorporate the 
new closure for the highly compressible turbulent cascade into 
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fluid mechani cs with adaptively refined large eddy simulations 
(FEARLESS: iMaier et alJl2009h . 
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